Theory of pattern-formation of metallic microparticles in poorly conducting liquid 
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We develop continuum theory of self-assembly and pattern formation in metallic microparticles 
immersed in a poorly conducting liquid in DC electric field. The theory is formulated in terms of two 
conservation laws for the densities of immobile particles (precipitate) and bouncing particles (gas) 
coupled to the Navier- Stokes equation for the liquid. This theory successfully reproduces correct 
topology of the phase diagram and primary patterns observed in the experiment [Sapozhnikov et al, 
Phys. Rev. Lett. 90, 114301 (2003)]: static crystals and honeycombs and dynamic pulsating rings 
and rotating mult i- petal vortices. 
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Electromagnetic manipulation and assembly of small 
particles in electrolytes has constitutes one of the great 
hopes for nanotechnology 1] and new generation of mi- 
cronuidic devices |2|. Many industrial technologies face 
the challenge of handling such single- or multi-component 
micro and nano- size ensembles. The dynamics of con- 
ducting microparticles in electric field in the air was stud- 
ied in 0, 3 • Phase transitions and clustering instability 
of the electrostatically driven conducting microparticles 
were found. In our recent work we reported new dynamic 
phenomena occurring with microparticles in poorly con- 
ducting liquid subject to strong electric field (up to 20 
kV/cm) p. It was shown that small metallic particles 
immersed in a toluene-ethanol mixture in DC electric 
field form a rich variety of novel phases. These phases 
include static honeycombs and two-dimensional crystals; 
dynamic multi-petal vortices and pulsating rings. The 
phenomena were attributed to interaction between par- 
ticles and electro-hydrodynamic (EHD) flows. 

In this Letter we develop theory of pattern forma- 
tion and self-assembly of metallic microparticles in a 
poorly conducting liquid placed between two horizontal 
planar electrodes. The theory is formulated in terms of 
two depth- averaged conservation laws for the densities of 
immobile particles (precipitate) and bouncing particles 
(gas) coupled to the Navier-Stokes equation for vertical 
component of velocity of the liquid v z , whereas horizon- 
tal velocity v± is obtained from the continuity condi- 
tion. This theory successfully reproduces topology of the 
phase diagram and primary patterns observed in the ex- 
periment: static crystals and honeycombs and dynamic 
pulsating rings and vortices Q . In the framework of our 
theory we demonstrate that the rotation of clusters is 
the result of a symmetry-breaking instability leading to 
formation of the travelling wave at the cluster perimeter. 

Model Schematics of experimental setting is shown in 
inset of Figure [I] Our studies show that two major con- 
trol parameters are the potential difference AU which 
determines average electric field E = —AU/d, d is the 
spacing between electrodes, and the concentration c of 



the additive (e.g. ethanol) which characterizes the con- 
ductivity of the liquid. Following the analysis of Ref. 
|j , we describe the evolution of particulate by the num- 
ber density of precipitate p p (r,t) and bouncing particles 
(or gas) p g (r,i), where r = (x,y) are horizontal coor- 
dinates. All the quantities are averaged over the ver- 
tical coordinate z. Since the total number of particles 
N = J (p p + p g )dxdy is conserved, the evolution of the 
particulate is described by the conservation laws 

d t p p = VJ p + / , d t Pg = \7Jg - f (1) 

Here J p ^ g are the mass fluxes of precipitate and gas re- 
spectively and the function / describes gas/precipitate 
conversion which depends on p Pj9 , electric field E and 
local concentration c. The fluxes can be written as: 

3p,g = D p ,gVp p ,g + OL pi g(E)V±P Pi g(l ~ ft '(E) Q p , g) (2) 

where D PyQ are precipitate/gas diffusivities. The last 
term, describing particles advection by fluid, is remi- 
niscent of Richardson- Zaki relation for drag force fre- 
quently used in the engineering literature [(j . The factor 
(1 — /3(E) p Pi9 ) describes saturation of the flux at large 
particle densities p ~ Experiments show that at 

the onset of motion the maximum density in the pat- 
terns such as honeycombs is below sub-monolayer cov- 
erage whereas for large values of E the particles form 
multi-layered structures. Thus (3 should decrease with 
the increase of E. According to the Richardson- Zaki re- 
lation, the coefficients a g ^ p decrease with the increase of 
p Q)P (i.e. void fraction). In order to mimic this effect on 
qualitative level we assumed that a g ^ p decreases with the 
increase of E (due to the increase of maximum density 
p ~ 1/(3 with E). Since the gas is mostly concentrated 
near the upper electrode and the precipitate near the 
bottom, gas and precipitate are advected in opposite di- 
rections, i.e. the transport coefficients a Py9 have opposite 
signs. Since gas is more mobile than precipitate, we set 
D g ^> D p . In the limit D g — > oo and for a P:g = one 
recovers the description of Ref. |4| . 



2 



We assume that vertical vorticity of the liquid Q z = 
d x v y — d y v x is small in comparison with the in-plane vor- 
ticity. This assumption is justified by the experimental 
observation that toroidal vortices create no or very small 
horizontal rotation. From Q z — one obtains 



-V0 



(3) 



where <p is a "quasipotential". Substituting Eq. (JHJ) 
into the continuity equation Vv = one expresses the 
quasipotential through vertical velocity 



V 2 = d z v z 



(4) 



The vertical velocity v z can be obtained form the cor- 
responding Navier-Stokes equation 



Po(d t v z + vV^) = vV 2 v z - d z p + E z q 



(5) 



where po is the density of liquid (we set po = 1), v is the 
viscosity, p is the pressure, and q is the charge density. 
The last term describes electric force acting on charged 
liquid. In order to average Eqs. (@J),(0 over the thickness 
of the cell < z < d, we assume that v z is symmetric 
with respect to d/2 and v± is antisymmetric. Then after 
the averaging and taking into account that d z v z = at 
z = 0, d one obtains 



d t V = vV 2 V - (V - Ap + (E z q) 



(6) 



where V — d~ x J Q d v z dz, and Ap = p(d) — p(0). Term 
— (V accounts for a small dissipation due to friction be- 
tween liquid, particles, and the walls of container. Using 
the symmetry condition and integrating Eq. (HJ) over the 
lower half of the cell (0 < z < d/2) one obtains 



V 2 $ = a V 



(7) 



where the constant ao ~ 0(1) can be scaled away and <E> 
is averaged <j). The volume charge density q ~ — c is nega- 
tive and proportional to additive concentration c, vertical 
component of the electric field E z depends on E and the 
local density of particles, p v + p g . Since the increase in 
the amount of conducting particles decreases the effective 
spacing between the electrodes and, therefore, increases 
the apparent electric field, one obtains 



E z 



E 



1 - (P P +Pg)/s 



E + E(p p + p g )/s (8) 



where constant s « d/ro, ro is particle diameter. After 
applying the divergence operator to the Navier-Stokes 
equation one finds that the pressure itself is a functional 
of (Eq). In the most general form it can be written 
as Ap = J K(r — r')(E z q)dr' . The kernel has property 
f Kdr' = 1 since uniform force distribution does not cre- 
ate net flow of the liquid. The precise form of the kernel 




FIG. 1: Real part of A vs k for three regimes: honeycomb 
(solid line, E — —45,(3 — 2,k = —0.1), coalescence (point- 
dashed line E = — 55, /3 — 0.5, ft = 0.1), stable 2D crystal 
(dashed line, E = 50, /3 = 2, ft = —0.1). Other parameters: 
C 0.02, p 0.3, D p = l,a p = a g = -0.3, v = 2,d = 
1,/io = 0. Inset: Schematics of experimental setup (right) and 
distance R between t- vortices vs t (solid line) for parameters 
of Figure 21 (left). Dashed line shows fit R ~ y^o — t. 



is not available on our level of description. We used the 
following kernel expressed by its Fourier transform 

K(k) = J exp\ikr}K(r)dr = exp[-K(E)k 2 - k 4 d%] (9) 

This form of the kernel is justified by the spatial 
isotropy (dependence of k 2 only), normalization (K(0) = 
J K(r')dr' = 1) and locality conditions (rapid decays for 
k > I/do? where do is the characteristic length of the or- 
der of electrodes separation). The field dependent factor 
k(E) describes experimentally observed transitions from 
honeycombs (ft < 0) at small fields to coalescence and 
attraction of toroidal vortices (ft > 0) for larger fields. 
After combining Eqs. (O,© and © one obtains 



d t V = vV 2 V -CV-cE J K^r- r f )(p p + p g )dr' (10) 

where the Fourier transform K\ = 1 — exp[— n(E)k 2 — 
k A d$]. We also scaled away s. 

The function / has a different structure for low con- 
centrations (c « 1) and high concentrations. For c — > 
this function coincides with that derived in Ref. 



fi = (pp-p*)- 



/ Ppi 



if < p p < 
Pg( l ~ Pp)i if P* < Pp < 1 • 



(11) 



where constants ~ 0(1) and p*(E) are discussed in 
Ref. 4]. For higher values of c the experiment suggests 
that both p g and p p tend to some equilibrium value which 
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is determined by the field E. In the linear approximation 
it can be written as 

f 2 = cC 2 (p g -p(E,q)p p ) (12) 

Here C2 characterizes relaxation time towards the equi- 
librium. In steady-state fa enforces the relation between 
the densities p g /p p — > /i. In turn, a depends on the field 
and local charge distribution q. We set for simplicity 
q = const. However, convective flows will affect the ions 
density, which will shift the gas/precipitate equilibrium, 
see 0. This effect can be modelled by including the de- 
pendence of a on V. We used the following expression 
a = /io(^)(tanh(/iisign(i£)V) + 1). The particular form 
of this function appears to be not relevant. It incorpo- 
rates the observation that gas concentration is suppressed 
by rising flows for E < due to excess of negative ions 
and vise versa. Finally, we used hybrid form for / valid 
for arbitrary c 

f = h + exp[-c/co]/i (13) 

where cq is some "crossover" concentration. 

Stability of homogeneous precipitate. In the limit of 
small electric fields E we can set a — > 0, i.e. equilibrium 
gas density p g — > 0. Then the stability of homogeneous 
precipitate (i.e. homogeneous "Wigner" crystal state in 
experiment) p p = p = const can be readily performed 
because the equation for p g splits off and the analysis is 
reduced to Eq. (|TU|) and Eq. Q for p p . We will focus 
on the case of high concentration c ^> Co and neglect the 
last term in Eq. (|13|) (case of c = was considered in 
Ref. 4]). In the linear order for periodic perturbations 
V, p p ~ exp[At + ikx] Eqs. (QJ, (fTU|) yield 

XV = -(vk 2 + ()V + E(l-exp[-K(E)k 2 -k 4 4])p p 
Xp p = -D p k 2 p p -a p (3p(l-(3p)V (14) 

The growth rate X(k) depends on the sign of the field 
E, see Figure ffl For negative E and ft < in a cer- 
tain parameter range A is positive in a narrow band near 
the optimal wavenumber ko , which is the hallmark of so- 
called short-wave instability. In this regime our solution 
indicates formation of stable honeycomb lattice with the 
scale determined by ko. For the same parameters, but for 
E > the homogeneous state p p = p is stable (dashed 
line on Figure Q}. This observation is consistent with ex- 
perimental fact that field reversal transforms honeycomb 
to homogeneous precipitate and vise versa. With the in- 
crease of E and ft the left edge of the unstable band of 
X(k) crosses zero and long-wave instability ensues. This 
regime corresponds to the coalescence of clusters. 

Qualitative phase diagram is shown in Figure [2j The 
positions of transition lines are approximate because field 
dependence of ft, /io, is not available at the moment. To 
describe the observed phenomenology, we only assumed 
that ft, /icb P* monotonically increase with the increase 
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AU = -Ed (arbitrary units) 



FIG. 2: Qualitative phase diagram, AU = —Ed is applied 
voltage (plus on top plate), c is concentration of additive. 
Domain sc denotes static clusters, wc denotes stable homo- 
geneous precipitate (Wigner crystals), he -honeycombs, utv 
and dtv up/down toroidal vortices. 
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FIG. 3: Formation of honeycomb, shown snapshots of p p , 
for p 0.3, v = 2,E = -50, a 0.02, a = 0,Ci = 5,c = 
1, Co — 0.1, C* = 10, Dp = I, do = 1, a p — —a g — —0.6, ft — 
—0.1, /3 = 2, domain of integration 80x80 dimensionless units. 
Images from left to right: t=10; 400; 2000. Black corresponds 
to p p = 0, white to max(yo p ). 

of E, and the functions a g ^ p: /3 decreases (1/(3 is limiting 
density of the t- vortex which increases with the increase 
of E). The line A = for short- wavelength instability for 
E < depicts the transition to honeycombs. The tran- 
sition from honeycombs to pulsating rings is identified 
as a transition from shortwave to long-wave instability, 
which roughly coincides with k(E) = 0. This transition 
is associated with an overall increase of gas concentration 
(increase of p(E)) and decrease of (3 (increase of maxi- 
mum density of p p ). The transition from static clusters 
to dynamic structures and honeycombs occurs approxi- 
mately with the increase of concentration c at c = c$. For 
E > the transition from stable precipitate to "down t- 
vortices" is approximately given by k(E) = 0. Note that 
in this case there is no distinction between "immobile" 
and "Wigner crystal" state because our model does not 
take into account dry friction and adhesion of particles 
to the bottom. 

Numerical solution of Eqs. ([!()[). ([T]) was performed by 
quasi-spectral method based on Fast Fourier Transforma- 
tion (FFT). Typically 256x256 FFT harmonics in peri- 
odic boundary conditions were used. For c > Co, E < 
and ft < corresponding to the case of short- wavelength 
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FIG. 4: Evolution of pulsating rings, shown snapshots of p p , 
for p 0.18,1/ = 2,E = -70, a 0.02, /x 0.25, /xi 

0. 5, Ci = 5,c = l,c = 0.1, C* = 10, D p = l,D g = 5,d = 

1, a p — — a g — —0.15, k — 0.1, /3 — 0.5, domain of integration 
80 x 80 units. Images from left to right: t=450; 1120; 3000. 




FIG. 5: Rotating down t- vortices, shown snapshots of p p , 
for p 0.18, v = 2,E = 160, a 0.02, fi 19.5, /xi 

0. 5, Ci = 0.45, c = 1, co = 0.1, C* = 10, D p = l,D g = 5, d = 

1. a p = — a p — —0.15, ft = 0.1, f3 = 0.5, domain of integration 
150 x 150 units. Images from left to right: t=16; 320; 540. 

instability we observed spontaneous formation of hon- 
eycomb lattice from random initial conditions, see Fig- 
ure 01 For the same parameters but E > the homo- 
geneous state p p = const (Wigner crystal) was stable. 
With the increase of E and for n > (i.e. in the case of 
long-wavelength instability) we observed the transition 
from honeycombs/ Wigner crystals to the regime of clus- 
ter attraction and coalescence. In this case EHD flows 
accelerate coarsening process of t-vortices dramatically. 
The inter-vortex distance R in this regime vs time is 
shown in Figure ^ inset. The fitting gives approximately 
R ^ (to — t) 1 ! 2 . It is qualitatively consistent with ex- 
perimental data. This behavior can be anticipated from 
solution of Eq. Q. In 2D t-vortices generate localized 
distributions of vertical velocity V, and, therefore, loga- 
rithmic quasipotentials ~ log |r — r^|, is the vortex 
position. Thus, the asymptotic horizontal velocity will 
decay as V± ~ l/|r — leading to R ~ (to — t) 1 / 2 . 

In the course of coalescence the vortices grow and be- 
come unstable, producing pulsating rings for E < (Fig- 
ure EJ) and rotating objects for E > (Figure EJ). We 
believe that the instability in both cases is caused by the 
same mechanism: dynamic coupling between particles 
density and vertical flows described by the function /. 
However, for E < 0, i.e for utv, the vortices are primar- 



ily built of a low-mobile precipitate phase with a small 
amount of gas above them. The instability occurs in the 
bulk of the vortex, resulting in formation and breaking 
up of "gas bubbles" inside the precipitate. For E > 
the vortices are formed by a mobile gas phase with a 
small content of precipitate. The instability occurs at the 
vortex edge yielding counter-propagating clock/anticlock 
waves of shape deformation. Eventually only one wave 
survives, creating the effect of rotation. 

In our experiments we observed also multi-petal vor- 
tices exhibiting almost solid-state rotation. These multi- 
petal vortices obviously generate horizontal flow of the 
liquid which is neglected in the framework of our theory. 
We believe that vertical vorticity is in fact the conse- 
quence rather than the reason of rotation. It is likely that 
the vortex shape instability present in our model eventu- 
ally triggers rotation of surrounding fluid and generates 
vertical vorticity. Explicit incorporation of the vertical 
vorticity will lead to substantial technical difficulties and 
likely will not change the qualitative features. 

In conclusion, we develop phenomenological continuum 
theory of pattern formation and self-assembly of metallic 
microparticles immersed in a poorly conducting liquid. 
This theory reproduces primary patterns observed in the 
experiment, and leads to an interesting prediction of the 
relation of rotation with the vortex edge instability. The 
parameters of the model can be extracted from molecu- 
lar dynamics simulations and generalization of "leaky di- 
electric model" 7] and validated by experiments. We are 
grateful to B. Meerson, Y. Tolmachev and W.-K. Kwok 
for useful discussions. This research was supported by 
the US DOE, grant W-31-109-ENG-38. 
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